/********************************************************************************
This do file generates Table 2, 4, 5, and C.1 presented in the paper.
********************************************************************************/
qui set matsize 11000
qui set maxvar 32767

cd "C:\Users\a-yamazaki\Dropbox\CGEValidation\written\JAERE_final"
use Carbone_et_al.dta, clear

gen province = 0
replace province = 1 if GEO == "Alberta"
replace province = 2 if GEO == "British Columbia"
replace province = 3 if GEO == "Manitoba"
replace province = 4 if GEO == "New Brunswick"
replace province = 5 if GEO == "Newfoundland and Labrador"
replace province = 6 if GEO == "Northwest Territories" /*I dropped this from the estimation*/
replace province = 7 if GEO == "Nova Scotia"
replace province = 8 if GEO == "Nunavut"
replace province = 9 if GEO == "Ontario"
replace province = 10 if GEO == "Prince Edward Island"
replace province = 11 if GEO == "Quebec"  /*I dropped this from the estimation*/
replace province = 12 if GEO == "Saskatchewan"
replace province = 13 if GEO == "Yukon"

/* Variable labels: 
intensity07 is the EI_i tilde in EI_{ijr} equation 
CGEemission# are the EI_{jr} bar in EI_{ijr} equation. 
There are four versions as we explore different definitions of emission intensity

CGEemission --> direct+indirect emission/Gross Output
CGEemission2 --> direct+indirect/GDP
CGEemission3 --> direct emission/Gross Output
CGEemission4 --> direct emission/GDP
*/

gen CGEemission = 0
replace CGEemission =	0.144	if CGE == "AGR" & CGEemission !=.
replace CGEemission =	1.193	if CGE == "CEM" & CGEemission !=.
replace CGEemission =	0.446	if CGE == "CHM" & CGEemission !=.
replace CGEemission =	1.051	if CGE == "COL" & CGEemission !=.
replace CGEemission =	0.039	if CGE == "CON" & CGEemission !=.
replace CGEemission =	0.775	if CGE == "CRU" & CGEemission !=.
replace CGEemission =	0.169	if CGE == "ELE" & CGEemission !=.
replace CGEemission =	0.497	if CGE == "GAS" & CGEemission !=.
replace CGEemission =	0.046	if CGE == "GOV" & CGEemission !=.
replace CGEemission =	0.079	if CGE == "MFR" & CGEemission !=.
replace CGEemission =	0.272	if CGE == "MIN" & CGEemission !=.
replace CGEemission =	0.368	if CGE == "OIL" & CGEemission !=.
replace CGEemission =	0.329	if CGE == "PPP" & CGEemission !=.
replace CGEemission =	0.619	if CGE == "PRM" & CGEemission !=.
replace CGEemission =	0.042	if CGE == "SER" & CGEemission !=.
replace CGEemission =	0.051	if CGE == "TRD" & CGEemission !=.
replace CGEemission =	0.666	if CGE == "TRN" & CGEemission !=.

cap gen CGEemission2 = 0
replace CGEemission2 =	0.573	if CGE == "AGR" & CGEemission2 !=.
replace CGEemission2 =	8.845	if CGE == "CEM" & CGEemission2 !=.
replace CGEemission2 =	1.049	if CGE == "CHM" & CGEemission2 !=.
replace CGEemission2 =	1.559	if CGE == "COL" & CGEemission2 !=.
replace CGEemission2 =	0.323	if CGE == "CON" & CGEemission2 !=.
replace CGEemission2 =	5.621	if CGE == "CRU" & CGEemission2 !=.
replace CGEemission2 =	0.317	if CGE == "ELE" & CGEemission2 !=.
replace CGEemission2 =	0.769	if CGE == "GAS" & CGEemission2 !=.
replace CGEemission2 =	0.166	if CGE == "GOV" & CGEemission2 !=.
replace CGEemission2 =	0.675	if CGE == "MFR" & CGEemission2 !=.
replace CGEemission2 =	0.979	if CGE == "MIN" & CGEemission2 !=.
replace CGEemission2 =	7.516	if CGE == "OIL" & CGEemission2 !=.
replace CGEemission2 =	0.997	if CGE == "PPP" & CGEemission2 !=.
replace CGEemission2 =	2.34	if CGE == "PRM" & CGEemission2 !=.
replace CGEemission2 =	0.165	if CGE == "SER" & CGEemission2 !=.
replace CGEemission2 =	0.169	if CGE == "TRD" & CGEemission2 !=.
replace CGEemission2 =	3.244	if CGE == "TRN" & CGEemission2 !=.

cap gen CGEemission3 = 0
replace CGEemission3 =	0.097	if CGE == "AGR" & CGEemission3 !=.
replace CGEemission3 =	1.054	if CGE == "CEM" & CGEemission3 !=.
replace CGEemission3 =	0.396	if CGE == "CHM" & CGEemission3 !=.
replace CGEemission3 =	0.996	if CGE == "COL" & CGEemission3 !=.
replace CGEemission3 =	0.009	if CGE == "CON" & CGEemission3 !=.
replace CGEemission3 =	0.681	if CGE == "CRU" & CGEemission3 !=.
replace CGEemission3 =	0.149	if CGE == "ELE" & CGEemission3 !=.
replace CGEemission3 =	0.461	if CGE == "GAS" & CGEemission3 !=.
replace CGEemission3 =	0.017	if CGE == "GOV" & CGEemission3 !=.
replace CGEemission3 =	0.034	if CGE == "MFR" & CGEemission3 !=.
replace CGEemission3 =	0.22	if CGE == "MIN" & CGEemission3 !=.
replace CGEemission3 =	0.285	if CGE == "OIL" & CGEemission3 !=.
replace CGEemission3 =	0.267	if CGE == "PPP" & CGEemission3 !=.
replace CGEemission3 =	0.541	if CGE == "PRM" & CGEemission3 !=.
replace CGEemission3 =	0.011	if CGE == "SER" & CGEemission3 !=.
replace CGEemission3 =	0.018	if CGE == "TRD" & CGEemission3 !=.
replace CGEemission3 =	0.381	if CGE == "TRN" & CGEemission3 !=.

cap gen CGEemission4 = 0
replace CGEemission4 =	0.229	if CGE == "AGR" & CGEemission4 !=.
replace CGEemission4 =	3.584	if CGE == "CEM" & CGEemission4 !=.
replace CGEemission4 =	0.796	if CGE == "CHM" & CGEemission4 !=.
replace CGEemission4 =	1.241	if CGE == "COL" & CGEemission4 !=.
replace CGEemission4 =	0.022	if CGE == "CON" & CGEemission4 !=.
replace CGEemission4 =	3.538	if CGE == "CRU" & CGEemission4 !=.
replace CGEemission4 =	0.208	if CGE == "ELE" & CGEemission4 !=.
replace CGEemission4 =	0.566	if CGE == "GAS" & CGEemission4 !=.
replace CGEemission4 =	0.03	if CGE == "GOV" & CGEemission4 !=.
replace CGEemission4 =	0.127	if CGE == "MFR" & CGEemission4 !=.
replace CGEemission4 =	0.394	if CGE == "MIN" & CGEemission4 !=.
replace CGEemission4 =	3.209	if CGE == "OIL" & CGEemission4 !=.
replace CGEemission4 =	0.605	if CGE == "PPP" & CGEemission4 !=.
replace CGEemission4 =	1.384	if CGE == "PRM" & CGEemission4 !=.
replace CGEemission4 =	0.018	if CGE == "SER" & CGEemission4 !=.
replace CGEemission4 =	0.029	if CGE == "TRD" & CGEemission4 !=.
replace CGEemission4 =	1.131	if CGE == "TRN" & CGEemission4 !=.

/*** Control variables 
Oil price
World trade index
US unemployment 
population by ages 
**************************************/
tab GEO, gen(province_)
sort year
by year:gen trends = 1 if _n == 1
replace trends = sum(trends)

forvalue i = 1(1) 11{
gen provincetrend`i' = province_`i'*trends
gen provincesqtrend`i' = provincetrend`i'*provincetrend`i'
}

*** eia_oil_price.dta ***** 
label var oil_price "EIA oil price"
gen lnOil = ln(oil_price)
forvalue i = 1(1)11{
gen oil_GEO_`i' = lnOil*province_`i'
}
label var US_unemp "US unemployment rate"
** from Labor Force Statistics from the Current Population Survey (LNS14000000)
forvalue i = 1(1)11{
gen USunemp_GEO_`i' = US_unemp*province_`i'
}

*** Export volume relative to 1913 **** 
** https://ourworldindata.org/trade-and-globalization ** 
gen WX = 0
replace WX = 19.24 if year == 2001
replace WX = 19.33 if year == 2002
replace WX = 20.19 if year == 2003
replace WX = 21.88 if year == 2004
replace WX = 22.93 if year == 2005
replace WX = 24.27 if year == 2006
replace WX = 24.87 if year == 2007
replace WX = 26.23 if year == 2008
replace WX = 21.45 if year == 2009
replace WX = 21.45 if year == 2010
replace WX = 25.12 if year == 2011
replace WX = 24.76 if year == 2012
replace WX = 24.51 if year == 2013
label var WX "World export share"

forvalue i = 1(1)11{
gen WX_GEO_`i' =  WX*province_`i'
}

egen ID_province = group(province SECTOR) /* Sector here is industries */
egen CGE_province = group(province CGE)  /* CGE here is sector */


xtset ID_province year

** population data is from pop.dta (Table 1710000501) **
global agestructure "lnkids lnteen2 lntwenty lnthirty lnfourty lnfifty lnabove60"

****** Table 2 *************************************************************
qui xtreg lnL CGEieffect BCtax i.ID#i.year, fe cluster(CGE_province) 
est sto one 

qui xtreg lnL CGEieffect BCtax i.ID#i.year oil_GEO_* USunemp_GEO_* WX_GEO_* $agestructure, fe cluster(CGE_province)
est sto two

qui xtreg lnL CGEieffect BCtax i.ID#i.year provincetrend* provincesqtrend* oil_GEO_* USunemp_GEO_* WX_GEO_* $agestructure, fe cluster(CGE_province)
est sto three

qui xtreg lnL CGEieffect i.ID#i.year i.province#i.year, fe cluster(CGE_province) 
est sto four 

esttab one two three four, keep(CGEieffect BCtax) b(%9.4f) star(* 0.1 ** 0.05 *** 0.01) r2 se noconstant 

****** Table 5 *************************************************************
qui xtreg lnL CGEieffect TIeffect BCtax i.ID#i.year provincetrend* provincesqtrend* oil_GEO_* USunemp_GEO_* WX_GEO_* $agestructure, fe cluster(CGE_province)
est sto trade 

esttab three trade, keep(CGEieffect TIeffect BCtax) order(CGEieffect TIeffect BCtax) b(%9.4f) star(* 0.1 ** 0.05 *** 0.01) r2 se noconstant 

****** Table C.1 *************************************************************
*** Table 5 in v9 ******************************************************************************************************* 
global EI "CGEemission CGEemission2 CGEemission3 CGEemission4"

gen idenomenator = L*intensity07 
bys GEO year CGE: egen denomenator = sum(idenomenator)
cap bys GEO year CGE: egen CGE_L = sum(L) if L !=.
label var CGE_L "L at CGE aggregation"
foreach var of global EI{
gen numerator`var' = CGE_L*`var'
gen weight`var' = numerator`var'/denomenator
gen intensity`var' = intensity07*weight`var'
}
	
foreach var of global EI{
replace CGEieffect = intensity`var'*BCtax

qui xtreg lnL CGEieffect BCtax i.ID#i.year provincetrend* provincesqtrend* oil_GEO_* USunemp_GEO_* WX_GEO_* lnkids lnteen2 lntwenty lnthirty lnfourty lnfifty lnabove60, fe cluster(CGE_province)
est sto one_`var' 
}

esttab one_CGEemission one_CGEemission2 one_CGEemission3 one_CGEemission4, ///
keep(CGEieffect BCtax) star(* 0.1 ** 0.05 *** 0.01) b(%9.4f) r2 se noconstant 


****** Table 4 *************************************************************
insheet using "trade_intensity.csv", clear
save "trade_intensity.dta", replace

insheet using ghg_intensity.csv, clear 
save "pseudo_intensity.dta", replace

insheet using "all_result_emp.csv", clear

rename v1 scenario				
rename v2 iden
rename v3 sector
rename v4 region
rename v5 value	

replace sector = upper(sector)
replace region = lower(region)
replace iden = lower(iden)
gen ctax = 0
replace ctax = 10 if scenario == "bc10" & region == "bc"
replace ctax = 15 if scenario == "bc15" & region == "bc"
replace ctax = 20 if scenario == "bc20" & region == "bc"
replace ctax = 25 if scenario == "bc25" & region == "bc"
replace ctax = 30 if scenario == "bc30" & region == "bc"

gen year = 2007
replace year = 2008 if scenario == "bc10"
replace year = 2009 if scenario == "bc15"
replace year = 2010 if scenario == "bc20"
replace year = 2011 if scenario == "bc25"
replace year = 2012 if scenario == "bc30"

merge m:1 sector using "pseudo_intensity.dta"
drop _merge

gen ly = log(value)
gen intact = ghgint * ctax

drop if iden !="labour"
drop if region == "all"
egen xtid = group(sector region)
xtset xtid year

label var ctax "Tax"
label var intact "Tax $\times$ intensity"

egen sectorID = group(sector)

qui xtreg ly ctax intact i.year#i.sectorID, fe vce(cluster xtid)
est sto one4
	
* Now replace all non-bc values with bc value
gen newvalue = ly if scenario == "bmk" & region != "bc"
egen totalvalue = total(newvalue), by(sector region)
replace totalvalue = ly if region == "bc"

qui xtreg totalvalue ctax intact i.year#i.sectorID, fe vce(cluster xtid)
est sto two4

merge m:1 sector using "trade_intensity.dta"
gen tradeintact = tradeint * ctax
label var tradeintact "Tax $\times$ trade intensity"

qui xtreg totalvalue ctax intact tradeintact i.year#i.sectorID, fe vce(cluster xtid)
est sto three4

esttab one4 two4 three4, keep(intact ctax tradeintact) order(intact tradeintact ctax) b(%9.4f) star(* 0.1 ** 0.05 *** 0.01) r2 se noconstant 








